Universality classes in the random-storage sandpile model 
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The avalanche statistics in a stochastic sandpile model where toppling takes place with a probabil- 
ity p is investigated. The limiting case p — 1 corresponds to the Bak-Tang-Wiesenfeld (BTW) model 
with deterministic toppling rule. Based on the moment analysis of the distribution of avalanche sizes 
we conclude that for < p < p c the model belongs to the DP universality class while for p c < p < 1 
it belongs to the BTW universality class, where p c is identified with the critical probability for 
directed percolation in the corresponding lattice. 
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Sandpile automata were proposed as a paradigm of 
self-organized critical (SOC) phenomena (D. These sim- 
ple models capture its essential dynamics, which takes 
place in the form of avalanches of all sizes. At the early 
state of SOC theory it was believed that the critical state 
of sandpile automata is insensitive to changes in model 
parameters, however some recent works contradict this 
statement. For instance, Vespignani and Zapperi have 
shown that driving and dissipation rates actually act as 
control parameters, criticality is obtained after fine tun- 
ing of these fields. On the other hand, we have recently 
shown that a class of models with stochastic rules dis- 
play a transition from SOC to directed percolation (DP) 
with increasing the degree of stochasticity ||. Never- 
theless, before we make our final conclusion, we have to 
investigate if the original Bak-Tang-Wiesenfeld (BTW) 
sandpile automaton and these modified models belong 
to the same universality class, otherwise they would just 
be different models. One may thus ask: do determin- 
istic and stochastic sandpile models belong to the same 
universality class? 

The BTW automaton is defined on d-dimensional hy- 
percubic lattice of linear size L. On each site i an in- 
teger variable Zi is defined, which we call energy fol- 
lowing 0. Energy is added to the system by selecting 
a site at random and increasing its energy by one, i.e. 
Zi — > Zi + 1. The energy addition continues until one of 
the sites reaches or exceeds a threshold z c = 2d, then this 
site topples transferring energy to its neighbors. Toppling 
is defined by the set of rules Zi — > Zi — z c and Zj — ► Zj + 1 
at each of the j nearest neighbors. The toppling at one 
site may induce an avalanche of toppling events. Af- 
ter the avalanche has stopped we re-start adding energy 
to the system. Open boundary conditions are assumed. 
The toppling rules of this model are deterministic and 
the only randomness is introduced through the random 
energy addition. 

May be the simplest stochastic sandpile model is the 
Manna or d-state model In this case z c — d < 2d 
and, therefore, only d of the 2d neighbors will receive 
energy after toppling. These d neighbors are selected at 
random. Real-space renormalization group approach ||] 



suggests that the BTW and Manna models belongs to 
the same universality class. However, there is not com- 
plete agreement in numerical simulations. Early large 
scale simulations of the Manna and BTW models || 
in two dimensions, show that the avalanche distributions 
are described by the same exponents for the power law 
decay and the scaling of the cutoffs. These results where 
supported by more recent estimates of the avalanche ex- 
ponents by Liibeck and Usadel . On the contrary, Ben 
Hur and Biham || analyzed the scaling of conditional 
expectation values of various quantities, obtaining sig- 
nificant differences in the exponents for the two models. 
However, Chessa et al have recently shown || that the 
method of conditional expectation values, introduced by 
Ben Hur and Biham, is systematically biased by non- 
universal corrections and, therefore, does not provide in- 
dications on universality classes. Moreover, according to 
the large scale simulations performed by Chessa at al || , 
the BTW and Manna models belong to the same univer- 
sality class. 

In the Manna model the randomness appears in the 
energy transfer after toppling, but the condition for top- 
pling is deterministic. Motivated on a directed model by 
Tadic and Dhar [jl0| we have proposed a different model, 
where sites topples when z > z c = 2d but with proba- 
bility p Q. In this case, as in the BTW model, each of 
the 2d neighbors receives one energy grain but the con- 
dition for toppling is stochastic. In the particular case 
p = 1 one recovers the BTW model. For p < 1 sites 
accumulate a random amount of energy, random storage 
model (RSM) jO]. From mean- field theory and numer- 
ical simulations in one dimension || we obtain that the 
RSM self-organizes into a critical state for p c < p < 1, 
while it is subcritical for p < p c . The subcritical state is 
found to be similar to DP and p c is identified with the 
critical threshold for DP. The correlation length expo- 
nents are identical to those of DP but other exponents 
result different due to open boundary conditions. How- 
ever, in that occasion our analysis was limited to a small 
region close to p c and simulations were performed only in 
d = 1. Thus, we could not provide any indication about 
the universality classes. 
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In the present work we extend the numerical simula- 
tions of the RSM to the whole range of p, in one and two 
dimensions. To determine the avalanche size and dura- 
tion exponents we use moment analysis, a technique pre- 
viously introduced by de Menech at al |l^] to obtain the 
scaling exponents for the BTW model, techniques. The 
numerical evidence suggests that the RSM and the BTW 
model belong to different universality classes. Moreover, 
for p < p c we corroborate that the model is similar to 
DP. 

The MF theory of the RSM (| reveals that there is a 
critical probability p c below which the system is subcrit- 
ical. In this regime p is small and sites accumulates a 
large amount of grains. Let us assume that all sites have 
heights above the critical threshold z c . Then when a new 
grain is added at certain site this site will topples with 
probability p. If it topples then each of its nearest neigh- 
bors will also topple with probability p and so on. This 
picture is equivalent to a site directed percolation prob- 
lem, where sites become active with probability p if one 
of its neighbors was active in the previous step. Hence 
the correlation length is given by £ ~ (p c — p)~ v , where 
p c is the critical threshold for site directed percolation 
in the corresponding lattice and v is the spatial correla- 
tion length exponent. Now, if p <C p c in such a way that 
£ -C L, where L is the system size, then an avalanche 
starting from a bulk site will never reach the boundary 
and, therefore, no grain will be dissipated. Hence, the av- 
erage height will increase in time because dissipation will 
no balance the grain addition from the external field. Af- 
ter a time long enough most sites will have heights above 
z c , supporting our starting assumption. The dynamical 
state below p c is thus equivalent to DP and p c is identified 
with the site DP threshold in the corresponding lattice. 
Above p c and for z > z c at all sites the system will still 
be equivalent (at least for small times) to DP and an in- 
finite avalanche will be generated. However, such state 
is not stable because the infinite avalanche will reach the 
boundary leading to dissipation of grains and hence the 
decrease of the column heights up to the avalanche has 
stopped. Above p c the stationary state is no more equiv- 
alent to DP. The balance between the grains added by 
the external field and dissipation at the boundary leads 
to a SOC state. Our task will be to determine if this 
SOC state belongs to the same universality class of the 
BTW model. 

Let s be the avalanche size and T its duration, s is 
the number of toppling events in an avalanche. T is the 
number of steps required to obtain a stable configuration 
starting from the initial site with z > z c , which triggered 
the avalanche, taking into account that on each step all 
sites are updated in parallel following the toppling rule. 
Both, s and T, are random variables. Their distributions 
are given by, assuming scaling, 

P x {x) =x- r *f x (x/x c ), (1) 

where x is s or T, s c and T c are the cutoff avalanche 



size and duration, and f x is a cutoff function. The cut- 
offs scale with the characteristic length of the system £ 
according to 

z c ~£ & , (2) 

where (3 S = D, [3 t = z, D is the avalanche dimension 
and z the dynamic scaling exponent. In the SOC state 
£ ~ L and, therefore, the cutoffs scale with system size. 
On the contrary, in the subcritical state £ ~ (p c — p)~ u 
and, hence, the cutoffs diverges when p approaches p c . 
The scaling exponents t x and (3 X are not all independent. 
From the identity P(s)ds = P(T)dt one obtains 

(t. - 1)D = {n - l)z. (3) 

On the other hand, in the SOC state (p > p c ) the mean 
avalanche size should scale as (s) ~ L? and, therefore, 

(2 - t s )D = 2. (4) 

These scaling relations are useful to test the reliability of 
the numerical estimates. 

The purpose of present numerical simulations is to de- 
termine the exponents r s , r t , D and z. We are going 
to take up this task using as fundamental technique the 
moment analysis Jl2[ . The q moment is given by 

(x q ) = J dxP(x)x q (5) 

where 

<Tx{q) = Ac(l - t x ) + /3 x q. (6) 

The last equivalence in eq. (|^) is not valid for small 
values of q. For small q the integral depends in the func- 
tional form of P{x) in the whole range of x, while scaling 
assumptions are in general not valid for small x. The 
extreme case is q = 0, normalization imposes c^O) = 
and, therefore, eq. (||) is not valid. But, for q > t x — 1 
large x dominates leading to the last equivalence in eq. 
©■ 

After computing the moments one can obtain a x (q) 
from a liner fit to the log-log plot of (x q ) vs. £. Then one 
can obtain t x and (3 X from a linear fit to the straight part 
of the plot <J x (q) vs. q. Above p c we have £ ~ L, but it is 
a function of p below. To compute the correlation lenght 
below p c we use the following expression M 

Eoc ( • ■ \1 

t=0 l^i=Q\ l ~ l 0) Pai (n\ 

where io is the position of the initial active (z > z c ) site, 
t is the number of steps measured in the time scale of the 
avalanche and p a i — 1 {p a i = 0) in active (inactive) sites. 
Moreover, from the log-log plot of £ vs. p c — p one can 
obtain a numerical estimate of p c and v. 

d = 1: The BTW model (p — 1) in one dimension ex- 
hibits trivial critical behavior jjj . No power law behavior 
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is observed in the avalanche size and duration distribu- 
tions, however one can compute the exponents D and z 
from the scaling of the moments with system size, result- 
ing D 2 and z « 1 . Using this values and the scaling re- 
lations in eqs (J3J) and (Q) one could obtain the exponents 
t s and r t (t s = r t = 1) however these scaling laws are not 
valid in this case because the distributions of avalanche 
size and duration do not follow the scaling law in eq. 
(|l|). On the contrary, the RSM in one dimension has 
non trivial behavior. The q dependence of a s and at, for 
different values of p < 1, is shown in figs. |l|and|^, respec- 
tively. The numerical estimates of the scaling exponents 
are given in table |. For p c < p < 0.8 we observe that 
a s (q), D, t s and r t are practically insensitive to changes 
in p, but systematic deviations are observed for at (q) and 
z. For p — 0.9 the scaling exponents are between those for 
p = 1 and p = 0.8. On the other hand, using the numer- 
ical data below p c , we have obtained p c = 0.707 ± 0.002, 
v = 1.07 ± 0.03 and z = 1.57 ± 0.02 which are, within 
the numerical error, identical to the series expansion es- 
timates for site DP in a square lattice |Q. Moreover, 
the exponents S — r t — 1 and z are consistent with pre- 
vious numerical simulations ||. We have also carry out 
finite size scaling of the distributions of avalanche size 
and duration. In all cases, including the BTW limit, we 
observe a good data collapse and the obtained scaling ex- 
ponents are in agreement with those obtained from the 
moment analysis. Moreover, within the numerical error, 
the scaling relations in eqs. (||) and are satisfied. 

d = 2: The BTW model has nontrivial exponents. 
However, the data collapse was not compatible with the 
scaling assumption in eq. (|l|) . The corrections to scaling 
in the BTW has been found to be very strong , mak- 
ing necessary the use of very large lattice sizes to obtain 
accurate estimates of the scaling exponents. The largest 
lattice size used in our simulations, L = 512, seems to 
be not large enough. Using lattice sizes ranging from 
L = 512toL = 2048 Chessa et al |j] have obtained a 
good finite size scaling for the distribution of avalanche 
size, but have not for the distribution of avalanche du- 
ration. We thus rule out the possibility of determine 
the BTW exponents in two dimensions with such small 
lattice sizes. Instead of that, we are going to use the 
numerical estimates by Chessa et al ^ and Lubeck and 
Usadel |t| . On the contrary, the RSM displays good finite 
size scaling for the lattice sizes we have used. Moreover, 
the scaling exponents obtained from the finite size scaling 
are in agreement with those obtained from the moment 
analysis. The q dependence of a s (q) and at(q) is shown 
in figs. |^ and |j, respectively. The numerical estimates 
of the scaling exponents for different values of p < 1 are 
given in tab. [Fj], toguether with the reports by Chessa 
et al H and Lubeck and Usadel @ for the BTW model. 
a s (q), D, t s and r t are practically insensitive to changes 
in p above p c , even considering the deterministic limit. 
On the contrary, at(q) and z suffer from systematic de- 
viations with changing p. In this case we must be more 
careful because the finite size effects have stronger in- 



fluence on the distribution of avalanche duration. For 
instance, for the largest lattice size used L = 512 the 
distribution of avalanche sizes cover about six decades 
while the distribution of avalanche durations cover less 
than five decades. To obtain more precise determination 
of the dynamic scaling exponents we must increase sys- 
tem size. In the mean time, the scaling exponents of the 
distribution of avalanche sizes indicates that the RSM in 
the range p c < p < 1 belongs to the same universality 
class of the BTW model, which correspond with the limit 
p = 1. On the other hand, below p c we have obtained 
p c = 0.344 ±0.001, v = 0.728 ±0.002 and z = 1.76 ±0.02 
which are, within the numerical error, identical to the 
numerical estimates for site DP in a BCC lattice [Q. In 
this case the difference between the exponents below and 
above p c is significant, showing that the RSM above and 
below p c belongs to different universality classes. 

The numerical simulations in d = 2 corroborate that 
the RSM is similar to DP below p c . As in d = 1, the expo- 
nents v and z and the critical probability p c are identical 
to the DP values in the corresponding lattice. On the 
other hand, the difference between the exponents D and 
z, obtained using the data below p c and those obtained 
above but close to p c , is far form being contained within 
the error bars. The model above and below p c belong 
to different universality classes. Below p c it is DP with 
open boundaries and random initial seed, while above p c 
we will continue using the term RSM. 

The numerical evidence obtained for the avalanche size 
distribution indicates that in d = 2 the RSM belong to 
the universality class of the BTW model. The scaling 
exponents r s , r t and D are practically independent of 
p in the SOC regime p c < p < 1, however, z shows a 
strong p dependence which may be attributed to finite 
size effects. On the contrary in d = 1 the distribution 
of avalanche sizes and duration for the BTW model dis- 
play trivial behavior while in the RSM they satisfy the 
scaling hypothesis. In larger dimensions d > 2 we expect 
that the SOC regime of the RSM belongs to the BTW 
universality class as in d = 2 . 

In summary the RSM has three different regimes. I: 
< p < p c similar to DP, II: p c < p < 1 where the top- 
pling rule are still stochastic but the system is in a SOC 
state and III: p = 1 (BTW) the toppling rules are deter- 
ministic. Based on the moment analysis of the distribu- 
tion of avalanche sizes we conclude that for < p < p c 
the model belongs to the DP universality class while for 
p c < p < 1 it belongs to the BTW universality class. 

We thanks A. Vespignani for bringing up our attention 
on two-dimensional simulations. This work was partially 
supported by the Alma Mater prize, given by the Uni- 
versity of Havana. 
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TABLE I. Scaling exponents in d = 1 for different values 
of p. 
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TABLE II. Scaling exponents in d = 2 for different val- 
ues of p. * These exponents may be affected by finite size 
corrections. 
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FIG. 1. Plot of (J 3 (q) in d = 1 for different values of p. 
Data above p c — 0.707(2) where obtained using lattice sizes 
L — 80, 160, 320 and 640. Data below p c where obtained 
using probabilities p = 0.670, 0.688, 0.696 and 0.7000. 
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FIG. 3. Plot of <J s (q) in d = 2 for different values of p. 
Data above p c — 0.344(1) where obtained using lattice sizes 
L = 64, 128, 256 and 512. Data below p c where obtained 
using probabilities p = 0.26, 0.30, 0.33 and 0.34. 
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FIG. 2. Plot of <Jt{q) in d = 1. The lattice sites and prob- 
abilities used are the same as described in the caption of fig. 
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FIG. 4. Plot of (Jt(q) in d = 2. The lattice sites and prob- 
abilities used are the same as described in the caption of fig. 
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